rm(list=ls())

## Load Libraries

library(stargazer)
library(lubridate)
library(gridExtra)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(xts)
library(tidyquant)
library(lubridate) # for the wday() and ymd() functions
library(lfe)
library(plm)

# load data 

# data 1: twitter data
data <- readRDS("ps_data.RDS")
# data 2: twitter academic composition data 
top50 <- read_csv("ps_top50.csv")

########################################################################################## 
########################################################################################## 
############################## MAIN MANUSCRIPT ###########################################
########################################################################################## 
########################################################################################## 




##############################
########## FIGURE 1 #########
##############################

g <- ggplot(top50, aes(as.factor(Gender))) + geom_bar(aes(fill = Rank)) + 
  scale_fill_manual(
    values =c("#BCE5E6", "#2382A1", "#2A5078")) + theme_minimal() + 
  xlab("Gender") + ylab("Count")    +  
  theme(legend.title = element_blank()) + 
  theme(axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=10),
        axis.text.x = element_text(color= "black", size=10),
        axis.title.y = element_text(color="black", size=14, face="bold"),
        plot.title = element_text(size = 14, face="bold")) +
  ggtitle("All Faculty at Top 50 Schools (N=1,747)")  +  coord_cartesian(ylim =c(0, 1200))+  
  theme(legend.position = "none")

g <- g + ggplot2::annotate("text", x = 1, y = 90,fontface=2, label = "Full (12.7%)") +
  ggplot2::annotate("text", x = 1, y = 310, fontface=2, label = "Associate (10.48%)") +
  ggplot2::annotate("text", x = 1, y = 480, fontface=2, label = "Assistant (9.73%)") +
  ggplot2::annotate("text", x = 2, y = 320, fontface=2, label = "Full (39.84%)") +
  ggplot2::annotate("text", x = 2, y = 830, fontface=2, label = "Associate (16.88%)") +
  ggplot2::annotate("text", x = 2, y = 1070, fontface=2, label = "Assistant (10.42%)") 
g


g1 <- ggplot(top50[top50$twitteruser==1,], aes(as.factor(Gender))) + geom_bar(aes(fill = Rank)) + 
  scale_fill_manual(
    values =c("#BCE5E6", "#2382A1", "#2A5078")) + theme_minimal() + 
  xlab("Gender") + ylab("")    +  
  theme(legend.title = element_blank()) + 
  theme(axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=10),
        axis.text.x = element_text(color= "black", size=10),
        axis.title.y = element_text(color="black", size=14, face="bold"),
        plot.title = element_text(size = 14, face="bold")) +
  ggtitle("Twitter User (N= 725 (41.5%))") + scale_y_continuous(breaks = c(0, 250, 500, 750, 1000, 1250)) + 
  coord_cartesian(ylim =c(0, 1200)) +  
  theme(legend.position = "none")

g1 <- g1 + ggplot2::annotate("text", x = 1, y = 40, label = "Full (11.86%)", fontface=2) +
  ggplot2::annotate("text", x = 1, y = 125, label = "Associate (11.89%)", fontface=2) +
  ggplot2::annotate("text", x = 1, y = 220, label = "Assistant (14.21%)", fontface=2) +
  ggplot2:: annotate("text", x = 2, y = 100, label = "Full (30.76%)", fontface=2) +
  ggplot2:: annotate("text", x = 2, y = 280, label = "Associate (16.97%)", fontface=2) +
  ggplot2::annotate("text", x = 2, y = 400, label = "Assistant (14.62%)", fontface=2) 

g1



jpeg("twitteruser.jpeg", units="in", width=10, height=6.5, res=300)

egg::ggarrange(g, g1, nrow = 1)

dev.off()


##############################
############ FIGURE 2 ##########
##############################


dailyaverage <- data %>% group_by(gender,date) %>% dplyr::summarise(total=n(),family=mean(family,na.rm=T),
                                                                    work=mean(work,na.rm=T), .groups = 'drop')


dailyaverage <- data %>% group_by(gender,date) %>% 
  dplyr::summarise(total=n(),family=mean(family,na.rm=T),
                   work=mean(work,na.rm=T),
                   prez=mean(prez, na.rm=T),
                   teaching=mean(teaching, na.rm=T),
                   blm = mean(blm, na.rm=T), .groups = 'drop')
dailyaverage$female <- dplyr::recode(dailyaverage$gender, '1' = 'Female', '0' = 'Male')
dailyaverage$female <- as.factor(dailyaverage$female)
dailyaverage$date <- as.Date(dailyaverage$date )


dailyaverage$teaching <- dailyaverage$teaching*100
dailyaverage$blm <- dailyaverage$blm*100
dailyaverage$prez <- dailyaverage$prez*100
dailyaverage$family <- dailyaverage$family*100
dailyaverage$work <- dailyaverage$work*100
dailyaverage$gender <- as.factor(dailyaverage$gender)


#pdat6 <- summarySE(dailyaverage, measurevar="work", groupvars=c("date","gender"))


pdat6 <- group.means <-ddply(dailyaverage,c("date","gender"),summarise,mean=mean(work))
pdat6$work <- pdat6$mean

pdat6$date <- as.Date(pdat6$date)
pdat6$female <- dplyr::recode(pdat6$gender, '1' = "Female", '0' = "Male")
pdat6$female <- as.factor(pdat6$female )
plot6<-  ggplot(pdat6, aes(x=date, y=work, colour=female, group=female),
                scale_fill_manual(breaks = c("Female", "Male"), values =c("#5c1a33", "#749dae"))) +
  geom_jitter( size=2, aes(colour=female), alpha=0.5) + xlab('') + ylab('Proportion of Work-Related Tweets') +
  scale_color_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) +
  geom_smooth(aes(x=date, y=work, fill=female),span=.2) + 
  scale_fill_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) + 
  scale_x_date(date_breaks = "months" , date_labels = "%b'%y") +
  geom_vline(xintercept = as.numeric(ymd("2020-03-13")), linetype="solid", 
             color = "#405E79", size=0.5) +
  geom_vline(xintercept = as.numeric(ymd("2020-01-21")), linetype="solid", 
             color = "#405E79", size=0.5) +
  theme_minimal()+
  theme(legend.position = "top",   legend.justification='left',
        legend.title = element_blank(),
        plot.title = element_text(color="black", size = 15, face = "bold"),
        axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=11, face="bold"),
        axis.text.x = element_text(color= "black", size=11, face="bold"),
        axis.title.y = element_text(color="black", size=14, face="bold")) 

plot6 <- plot6 + 
  ggplot2::annotate("text", x = ymd("2020-01-21"), y = 17, fontface=2, size=3, 
                    label = "First COVID-19\n case in the US") +
  ggplot2::annotate("text", x = ymd("2020-03-13"), y = 17, fontface=2, size=3, 
                    label = "Trump declared\n a national emergency") 
plot6


pdat7 <- group.means<-ddply(dailyaverage,c("date","gender"),summarise,mean=mean(family))
pdat7$family <- pdat6$mean


pdat7 <- summarySE(dailyaverage, measurevar="family", groupvars=c("date","gender"))
pdat7$date <- as.Date(pdat7$date)
pdat7$female <- dplyr::recode(pdat7$gender, '1' = "Female", '0' = "Male")
pdat7$female <- as.factor(pdat7$female )
plot7<-  ggplot(pdat7, aes(x=date, y=family, colour=female, group=female),
                scale_fill_manual(breaks = c("Female", "Male"), values =c("#5c1a33", "#749dae"))) +
  geom_jitter(size=2, aes(colour=female), alpha=0.5) + xlab('') + ylab('Proportion of Family-Related Tweets') +
  scale_color_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) +
  geom_smooth(aes(x=date, y=family, fill=female),span=.2) + 
  scale_fill_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) + 
  scale_x_date(date_breaks = "months" , date_labels = "%b'%y") +
  geom_vline(xintercept = as.numeric(ymd("2020-03-13")), linetype="solid", 
             color = "#405E79", size=0.5) +
  geom_vline(xintercept = as.numeric(ymd("2020-01-21")), linetype="solid", 
             color = "#405E79", size=0.5) +
  theme_minimal()+
  theme(legend.position = "top",   legend.justification='left',
        legend.title = element_blank(),
        plot.title = element_text(color="black", size = 15, face = "bold"),
        axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=11, face="bold"),
        axis.text.x = element_text(color= "black", size=11, face="bold"),
        axis.title.y = element_text(color="black", size=14, face="bold")) 

plot7 <- plot7 + 
  ggplot2::annotate("text", x = ymd("2020-01-21"), y = 12, fontface=2, size=3, 
                    label = "First COVID-19\n case in the US") +
  ggplot2::annotate("text", x = ymd("2020-03-13"), y = 12, fontface=2, size=3, 
                    label = "Trump declared\n a national emergency") 


plot7


jpeg("dailyaverage.jpeg", units="in", width=10, height=12, res=300)

egg::ggarrange(plot6, plot7, nrow = 2)

dev.off()

##############################
############# TABLE 1 ########
##############################


m1 <- lm(twitteruser ~ female*assistant + female*full , data=top50)

table <- capture.output({stargazer(m1, 
                                   dep.var.labels=c("Twitter User" ),
                                   covariate.labels = c("Female", "Assistant Prof", "Associate Prof", "Female X Assistant", "Female x Associate"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   no.space=TRUE,
                                   star.char = c("+", "*", "**", "***"), 
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

##############################
############# TABLE 2 ########
##############################


data_weekly <- data
data_weekly$date <- as.Date(data_weekly$date)


data_weekly$date  <- ymd(data_weekly$date)
saturdays <- data_weekly[wday(data_weekly$date) == 7, ] # filter for Saturdays
startDate <- min(saturdays$date) # select first Saturday
data_weekly$week <- floor(as.numeric(difftime(data_weekly$date, startDate, units = "weeks")))

data_weekly$week[data_weekly$date=='2020-03-13']
# treatment week is week 40 if using June 1, 2019 - June 1 2020 data

profp_w <- data_weekly %>% 
  group_by(screen_name,week, gender, rank) %>% 
  dplyr::summarise(total=n(),family=mean(family,na.rm=T), work=mean(work, na.rm=T), 
                   prez=mean(prez, na.rm=T), blm=mean(blm, na.rm=T), teaching=mean(teaching, na.rm=T),
                   .groups = 'drop')

profp_w_n <- data_weekly %>% 
  group_by(screen_name,week, gender, rank) %>% 
  dplyr::summarise(total=n(),familyP=mean(family,na.rm=T), workP=mean(work, na.rm=T),
                   familyN=sum(family),workN=sum(work), .groups = 'drop')




profp_w$pre_post <- 0 
profp_w$pre_post[profp_w$week >= 40] <- 1 

profp_w$family <- profp_w$family*100
profp_w$teaching <- profp_w$teaching*100
profp_w$work <- profp_w$work*100 
profp_w$blm <- profp_w$blm*100
profp_w$prez <- profp_w$prez*100 

junior_w <- profp_w %>% filter(rank==1)
associate_w <- profp_w %>% filter(rank==2)
full_w <- profp_w %>% filter(rank==3)


### Models 

m1 <-  felm(family ~ gender:pre_post + total | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

m2 <-  felm(work ~ gender:pre_post+ total | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

m3 <-  felm(family ~ gender:pre_post+ total | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)

m4 <-  felm(work ~ gender:pre_post + total| factor(screen_name) + factor(week) | 0 | 
              screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)

m5 <-  felm(family ~ gender:pre_post + total| factor(screen_name) + factor(week) | 0 | 
              screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)

m6 <-  felm(work ~ gender:pre_post + total| factor(screen_name) + factor(week) | 0 | 
              screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)

m7 <-  felm(family ~ gender:pre_post+ total | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)

m8 <-  felm(work ~ gender:pre_post+ total | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)

stargazer(m1,m2,m3,m4,m5,m6,m7,m8,omit.stat = 'all',type='text')

table <- capture.output({stargazer(m1,m2,m3,m4,m5,m6, m7,m8, 
                                   dep.var.labels=rep(c('Family Tweet','Work Tweet'),4),
                                   covariate.labels = c('Total', 'Female*Pandemic'),
                                   dep.var.caption = "",
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   star.cutoffs = c(0.1, .05,.01,.001),
                                   column.labels = c('All Faculty','Junior','Associate', 'Full'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   label = "main",
                                   title = "The Pandemic Effect on Family- and Work-Related Tweets",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)  

##############################
############# TABLE 3 ########
##############################


# note: we used following words to classify teaching-related words 
#teaching_words <- paste(c('assign', 'assigning', 'assignment','class','course',
  #                        'educator','enroll','essay','grade','grading',
   #                       'hybrid','instruction','instruct', 'intro','lecture','module','office hour',
    #                      'prepping','professor','quiz','rec letter' , 'recommendation letter', 
     #                     'recording','recording','semester',
      #                    'service','session','slides','student',
       #                   'syllabi','syllabus','syllabus','synchronous','taught','teacher',
        #                  'teaching','textbook','undergrad'),collapse="|")



m1 <-  felm(teaching ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

m3 <-  felm(teaching ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)

m5 <-  felm(teaching ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)

m7 <-  felm(teaching ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
              screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)


stargazer(m1,m3, m5, m7, omit.stat = 'all',type='text')

table <- capture.output({stargazer(m1, m3, m5, m7, 
                                   dep.var.labels=rep(c('Teaching Tweet'),4),
                                   covariate.labels = c('Female*Pandemic'),
                                   dep.var.caption = "",
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   star.cutoffs = c(0.1, .05,.01,.001),
                                   column.labels = c('All Faculty','Assistant','Associate', 'Full'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   label = "main",
                                   title = "The Pandemic Effect on Family- and Work-Related Tweets",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)  





####################################################################################
####################################################################################
######################## APPENDIX ##################################################
####################################################################################
####################################################################################


######################################################################################################
#####  TABLE C1 : Table 2 without controlling for total tweets 
######################################################################################################


mod1 <-  felm(family ~ gender:pre_post   | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

mod2 <-  felm(work ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

mod3 <-  felm(family ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)


mod4 <-  felm(work ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)


mod5 <-  felm(family ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)


mod6 <-  felm(work ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)

mod7 <-  felm(family ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)


mod8 <-  felm(work ~ gender:pre_post  | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(mod1,mod2,mod3,mod4,mod5,mod6, mod7,mod8, 
                                   dep.var.labels=rep(c('Family Tweet','Work Tweet'),4),
                                   covariate.labels = c('Female*Pandemic'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   column.labels = c('All Faculty','Junior','Associate', 'Full'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)  



######################################################################################################
#####  TABLE C2 : Table 2 with controls for total tweet interactions 
######################################################################################################

mod1 <-  felm(family ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total  | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

mod2 <-  felm(work ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total| factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w,   exactDOF=TRUE,psdef=FALSE)

mod3 <-  felm(family ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)


mod4 <-  felm(work ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w,   exactDOF=TRUE,psdef=FALSE)


mod5 <-  felm(family ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)


mod6 <-  felm(work ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w,   exactDOF=TRUE,psdef=FALSE)

mod7 <-  felm(family ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total| factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)


mod8 <-  felm(work ~ gender:pre_post + total+ gender:total + pre_post:total + gender:pre_post:total | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(mod1,mod2,mod3,mod4,mod5,mod6, mod7,mod8, 
                                   dep.var.labels=rep(c('Family Tweet','Work Tweet'),4),
                                   covariate.labels = c('Total', 'Female*Pandemic', 'Female*Total', 'Pandemic*Total', 'Female*Pandemic*Total'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   column.labels = c('All Faculty','Junior','Associate', 'Full'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)  



######################################################################################################
#####  TABLE C3 : Table 2 with tweet counts as the outcome 
######################################################################################################

profp_w_n <- data_weekly %>% 
  group_by(screen_name,week, gender, rank) %>% 
  dplyr::summarise(total=n(),familyP=mean(family,na.rm=T), workP=mean(work, na.rm=T),
                   familyN=sum(family),workN=sum(work), .groups = 'drop')

profp_w_n$pre_post <- 0 
profp_w_n$pre_post[profp_w_n$week >= 40] <- 1 

junior_w_n<- profp_w_n %>% filter(rank==1)
associate_w_n <- profp_w_n %>% filter(rank==2)
full_w_n <- profp_w_n %>% filter(rank==3)


mod1 <-  felm(familyN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w_n,   exactDOF=TRUE,psdef=FALSE)


mod2 <-  felm(workN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w_n,   exactDOF=TRUE,psdef=FALSE)

mod3 <-  felm(familyN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w_n,   exactDOF=TRUE,psdef=FALSE)


mod4 <-  felm(workN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w_n,   exactDOF=TRUE,psdef=FALSE)


mod5 <-  felm(familyN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w_n,   exactDOF=TRUE,psdef=FALSE)


mod6 <-  felm(workN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w_n,   exactDOF=TRUE,psdef=FALSE)

mod7 <-  felm(familyN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w_n,   exactDOF=TRUE,psdef=FALSE)


mod8 <-  felm(workN ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w_n,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(mod1,mod2,mod3,mod4,mod5,mod6, mod7,mod8, 
                                   dep.var.labels=rep(c('Family Tweet','Work Tweet'),4),
                                   covariate.labels = c('Female*Pandemic'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   column.labels = c('All Faculty','Junior','Associate', 'Full'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)  


######################################################################################################
#####  TABLE C4 : Table 2 with alternative treatment date (using the first week of March as the cutoff)
######################################################################################################


profp_w2 <- data_weekly %>% 
  group_by(screen_name,week, gender, rank) %>% 
  dplyr::summarise(total=n(),family=mean(family,na.rm=T), work=mean(work, na.rm=T))


profp_w2$pre_post <- 0 
profp_w2$pre_post[profp_w$week >= 39] <- 1 

profp_w2$work <- profp_w2$work*100 
profp_w2$family <- profp_w2$family*100 

junior_w2 <- profp_w2 %>% filter(rank==1)
associate_w2 <- profp_w2 %>% filter(rank==2)
full_w2 <- profp_w2 %>% filter(rank==3)

#install.packages("lfe")
#install.packages("plm")



mod1 <-  felm(family ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w2,   exactDOF=TRUE,psdef=FALSE)


mod2 <-  felm(work ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = profp_w2,   exactDOF=TRUE,psdef=FALSE)

mod3 <-  felm(family ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w2,   exactDOF=TRUE,psdef=FALSE)


mod4 <-  felm(work ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = junior_w2,   exactDOF=TRUE,psdef=FALSE)


mod5 <-  felm(family ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w2,   exactDOF=TRUE,psdef=FALSE)


mod6 <-  felm(work ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = associate_w2,   exactDOF=TRUE,psdef=FALSE)

mod7 <-  felm(family ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w2,   exactDOF=TRUE,psdef=FALSE)


mod8 <-  felm(work ~ gender:pre_post | factor(screen_name) + factor(week) | 0 | 
                screen_name, data = full_w2,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(mod1,mod2,mod3,mod4,mod5,mod6, mod7,mod8, 
                                   dep.var.labels=rep(c('Family Tweet','Work Tweet'),4),
                                   covariate.labels = c('Female*Pandemic'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   column.labels = c('All Faculty','Junior','Associate', 'Full'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)  


######################################################################################################
#####  FIGURE D1 : Pre-Treatment Gender Difference in Work- and Family-Related Tweets
######################################################################################################

profp_pt <- profp_w %>% filter(week <= 40)
profp_pt$week <- relevel(as.factor(profp_pt$week), ref='40')

junior_pt  <- profp_pt %>% filter(rank==1)
associate_pt  <- profp_pt %>% filter(rank==2)
full_pt <- profp_pt %>% filter(rank==3)

library(plm)


f<-plm(family~gender*(relevel(as.factor(week), ref='40')),data=profp_pt, 
       model = 'within', index = c('screen_name'))


w<-plm(work~gender*(relevel(as.factor(week), ref='40')),data=profp_pt, 
       model = 'within', index = c('screen_name'))

table <- capture.output({stargazer(f,w, 
                                   dep.var.labels=rep(c('Family Tweet','Work Tweet'),4),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   column.labels = c('All Faculty'),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes")),
                                   digits=2,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

f_df <- tidy(f) 
f_df <- f_df[41:80,]
w_df <- tidy(w)
w_df <- w_df[41:80,]

f_df$term <- c("Female x Week:-40",
               "Female x Week:-39",
               "Female x Week:-38", 
               "Female x Week:-37",
               "Female x Week:-36", 
               "Female x Week:-35",
               "Female x Week:-34", 
               "Female x Week:-33",
               "Female x Week:-32", 
               "Female x Week:-31", 
               "Female x Week:-30", 
               "Female x Week:-29", 
               "Female x Week:-28", 
               "Female x Week:-27", 
               "Female x Week:-26", 
               "Female x Week:-25",
               "Female x Week:-24", 
               "Female x Week:-23", 
               "Female x Week:-22", 
               "Female x Week:-21", 
               "Female x Week:-20", 
               "Female x Week:-19", 
               "Female x Week:-18", 
               "Female x Week:-17", 
               "Female x Week:-16", 
               "Female x Week:-15", 
               "Female x Week:-14", 
               "Female x Week:-13", 
               "Female x Week:-12", 
               "Female x Week:-11", 
               "Female x Week:-10", 
               "Female x Week:-9", 
               "Female x Week:-8", 
               "Female x Week:-7", 
               "Female x Week:-6", 
               "Female x Week:-5", 
               "Female x Week:-4", 
               "Female x Week:-3", 
               "Female x Week:-2", 
               "Female x Week:-1")

w_df$term <- c("Female x Week:-40",
               "Female x Week:-39",
               "Female x Week:-38", 
               "Female x Week:-37",
               "Female x Week:-36", 
               "Female x Week:-35",
               "Female x Week:-34", 
               "Female x Week:-33",
               "Female x Week:-32", 
               "Female x Week:-31", 
               "Female x Week:-30", 
               "Female x Week:-29", 
               "Female x Week:-28", 
               "Female x Week:-27", 
               "Female x Week:-26", 
               "Female x Week:-25",
               "Female x Week:-24", 
               "Female x Week:-23", 
               "Female x Week:-22", 
               "Female x Week:-21", 
               "Female x Week:-20", 
               "Female x Week:-19", 
               "Female x Week:-18", 
               "Female x Week:-17", 
               "Female x Week:-16", 
               "Female x Week:-15", 
               "Female x Week:-14", 
               "Female x Week:-13", 
               "Female x Week:-12", 
               "Female x Week:-11", 
               "Female x Week:-10", 
               "Female x Week:-9", 
               "Female x Week:-8", 
               "Female x Week:-7", 
               "Female x Week:-6", 
               "Female x Week:-5", 
               "Female x Week:-4", 
               "Female x Week:-3", 
               "Female x Week:-2", 
               "Female x Week:-1")

w_df <- w_df %>% mutate(model = "Work")
f_df <- f_df %>% mutate(model = "Family")

two_models <- rbind(w_df, f_df)
library(dotwhisker)
dwplot(two_models)  
dwplot(w_df)

cols <- c("Work-Related Tweets" = "#273253"  , "Family-Related Tweets" = "#ED5C4D")

d <- dwplot(two_models) +
  theme_minimal() + xlab("Coefficient Estimate") + ylab("") +
  theme(plot.title = element_text(face="bold"),
        legend.title = element_blank(),
        axis.title.x = element_text(color="black", size=12),
        axis.text.y = element_text(color= "black", size=11),
        axis.text.x = element_text(color= "black", size=11),
        axis.title.y = element_text(color="black", size=12)) +
  scale_colour_manual(name = c("Work", "Family"), 
                      values = c("#273253", "#078CA9")) 

d <- d + geom_vline(xintercept = 0, colour = "#72132C", linetype = 3, size=0.7)  + coord_flip() +
  theme(axis.text.x = element_text(angle=90), 
        legend.position = c(0.01, 0.85),
        legend.justification = c(0, 0))

d


######################################################################################################
#####  FIGURE E1: Daily Trends in % Teaching-Related Tweets
######################################################################################################

pdat6 <- summarySE(dailyaverage, measurevar="teaching", groupvars=c("date","gender"))
pdat6$date <- as.Date(pdat6$date)
pdat6$female <- recode(pdat6$gender, '1' = 'Female', '0' = 'Male')
pdat6$female <- as.factor(pdat6$female )
plot6<-  ggplot(pdat6, aes(x=date, y=teaching, colour=female, group=female),
                scale_fill_manual(breaks = c("Female", "Male"), values =c("#5c1a33", "#749dae"))) +
  geom_jitter( size=2, aes(colour=female), alpha=0.5) + xlab('') + ylab('Proportion of Teaching-Related Tweets') +
  scale_color_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) +
  geom_smooth(aes(x=date, y=teaching, fill=female),span=.2) + 
  scale_fill_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) + 
  scale_x_date(date_breaks = "months" , date_labels = "%b'%y") +
  geom_vline(xintercept = as.numeric(ymd("2020-03-13")), linetype="solid", 
             color = "#405E79", size=0.5) +
  geom_vline(xintercept = as.numeric(ymd("2020-01-21")), linetype="solid", 
             color = "#405E79", size=0.5) +
  theme_minimal()+
  theme(legend.position = "top",   legend.justification='left',
        legend.title = element_blank(),
        plot.title = element_text(color="black", size = 15, face = "bold"),
        axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=11, face="bold"),
        axis.text.x = element_text(color= "black", size=11, face="bold"),
        axis.title.y = element_text(color="black", size=14, face="bold")) 

plot6 <- plot6 + 
  ggplot2::annotate("text", x = ymd("2020-01-21"), y = 16, fontface=2, size=3, 
                    label = "First COVID-19\n case in the US") +
  ggplot2::annotate("text", x = ymd("2020-03-13"), y = 16, fontface=2, size=3, 
                    label = "Trump declared\n a national emergency") 
plot6


jpeg("teaching.jpeg", units="in", width=10, height=7, res=300)

plot6

dev.off()

######################################################################################################
#####  FIGURE E2: Daily Trends in % Trump and Biden-Related Tweets
######################################################################################################

pdat6 <- summarySE(dailyaverage, measurevar="prez", groupvars=c("date","gender"))
pdat6$date <- as.Date(pdat6$date)
pdat6$female <- recode(pdat6$gender, '1' = 'Female', '0' = 'Male')
pdat6$female <- as.factor(pdat6$female )
plot6<-  ggplot(pdat6, aes(x=date, y=prez, colour=female, group=female),
                scale_fill_manual(breaks = c("Female", "Male"), values =c("#5c1a33", "#749dae"))) +
  geom_jitter( size=2, aes(colour=female), alpha=0.5) + xlab('') + ylab('Proportion of Trump and Biden-Related Tweets') +
  scale_color_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) +
  geom_smooth(aes(x=date, y=prez, fill=female),span=.2) + 
  scale_fill_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) + 
  scale_x_date(date_breaks = "months" , date_labels = "%b'%y") +
  geom_vline(xintercept = as.numeric(ymd("2020-03-13")), linetype="solid", 
             color = "#405E79", size=0.5) +
  geom_vline(xintercept = as.numeric(ymd("2020-01-21")), linetype="solid", 
             color = "#405E79", size=0.5) +
  theme_minimal()+
  theme(legend.position = "top",   legend.justification='left',
        legend.title = element_blank(),
        plot.title = element_text(color="black", size = 15, face = "bold"),
        axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=11, face="bold"),
        axis.text.x = element_text(color= "black", size=11, face="bold"),
        axis.title.y = element_text(color="black", size=14, face="bold")) 

plot6 <- plot6 + 
  ggplot2::annotate("text", x = ymd("2020-01-21"), y = 16, fontface=2, size=3, 
                    label = "First COVID-19\n case in the US") +
  ggplot2::annotate("text", x = ymd("2020-03-13"), y = 16, fontface=2, size=3, 
                    label = "Trump declared\n a national emergency") 
plot6


jpeg("president.jpeg", units="in", width=10, height=7, res=300)

plot6

dev.off()


######################################################################################################
#####  TABLE E3: The Lockdown Effect on Trump- and Biden-Related Tweets 
######################################################################################################



prez_profp_w <- data_weekly %>% 
  filter(as.Date(date)>=as.Date('2019-06-01') & as.Date(date)<as.Date('2020-06-01')) %>%
  group_by(screen_name,gender,rank,pre_post, week) %>% 
  dplyr::summarise(total=n(),prez=mean(prez,na.rm=T))

prez_profp_w$prez <- prez_profp_w$prez *100 

#prez_profp_w$pre_post <- 0 
#prez_profp_w$pre_post[profp_w$week >= 40] <- 1 

prez_junior_w <- prez_profp_w %>% filter(rank==1)
prez_associate_w <- prez_profp_w %>% filter(rank==2)
prez_full_w <- prez_profp_w %>% filter(rank==3)


mod1 <- felm(prez ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
               screen_name, data = prez_profp_w,   exactDOF=TRUE,psdef=FALSE)
mod3 <- felm(prez ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
               screen_name, data = prez_junior_w,   exactDOF=TRUE,psdef=FALSE)
mod7 <- felm(prez ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
               screen_name, data = prez_associate_w,   exactDOF=TRUE,psdef=FALSE)
mod9 <- felm(prez ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
               screen_name, data = prez_full_w,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(mod1,mod3,mod7,mod9, 
                                   dep.var.labels= c('Trump- and Biden- Related Tweets'),
                                   covariate.labels = c('Female*Pandemic'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = 'week' , 
                                   column.labels = c('All Faculty','Assistant','Associate', 'Full'), column.separate = c(1,1,1,1),
                                   no.space=TRUE,
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   add.lines = list(c("Individual Fixed Effect?", "Yes", "Yes","Yes","Yes","Yes","Yes","Yes","Yes"),
                                                    c("Time Fixed Effect", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


######################################################################################################
#####  FIGURE E4: Daily Trends in % Black Lives Matter Tweets
######################################################################################################

## we used folloiwng BLM-related keywords 

blm_words <-  paste(c('blm', 'black lives matter', 'blacklivesmatter', 
                      'black lives', 'arbery', 'ahmaud arbery', 'ahmaud marquez arbery',  
                      'breonnataylor','breonna taylor', 'breonnataylor',  'floyd', 'george floyd', 
                      'georgefloyd', 'derek chauvin', 'chauvin', 'police', 'police brutality', 'icantbreathe',  'minneapolis',  
                      'myles cosgrove', 'cosgrove', 'joshua jaynes', 'jaynes', 'gregory mcmichael',
                      'travis mcmichael', 'mcmichael', 'william roddie bryan', 'bryan', 'louisville'), collapse='|')

#oldhandle$username <- oldhandle$oldhandle
#oldhandle$dummy <- 1 
#final3 <- merge(final2, oldhandle, by="username", all=T)
#final4 <- final3 %>% filter(dummy==1)
#final4$blm <- ifelse(grepl(blm_words,final4$tweet),1,0)
#data2 <- final4 %>% filter(as.Date(date)>=as.Date('2020-01-01') & as.Date(date)<as.Date('2020-08-01'))
#write.csv(data2, "ps_blmdata.csv")

data2 <- read.csv("ps_blmdata.csv")

dailyaverage2 <- data2 %>% group_by(gender,date) %>% dplyr::summarise(total=n(),
                                                                      family=mean(family,na.rm=T),work=mean(work,na.rm=T),
                                                                      blm=mean(blm, na.rm=T), prez=mean(prez, na.rm=T),
                                                                     )
dailyaverage2 <- dailyaverage2  %>% filter(is.na(gender)==F)
dailyaverage2$blm <-    dailyaverage2$blm*100


pdat8 <- summarySE(dailyaverage2, measurevar="blm", groupvars=c("date","gender"))
pdat8$date <- as.Date(pdat8$date)
pdat8$female <- recode(pdat8$gender, '1' = 'Female', '0' = 'Male')
pdat8$female <- as.factor(pdat8$female )
plot8<-  ggplot(pdat8, aes(x=date, y=blm, colour=female, group=female),
                scale_fill_manual(breaks = c("Female", "Male"), values =c("#5c1a33", "#749dae"))) +
  geom_jitter( size=2, aes(colour=female), alpha=0.5) + xlab('') + ylab('Proportion of BLM-Related Tweets') +
  scale_color_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) +
  geom_smooth(aes(x=date, y=blm, fill=female),span=.2) + 
  scale_fill_manual(values = c("Female" = "#5c1a33", "Male" = "#749dae")) + 
  scale_x_date(date_breaks = "months" , date_labels = "%b'%y") +
  geom_vline(xintercept = as.numeric(ymd("2020-03-13")), linetype="solid", 
             color = "#405E79", size=0.5) +
  geom_vline(xintercept = as.numeric(ymd("2020-05-25")), linetype="solid", 
             color = "#405E79", size=0.5) +
  theme_minimal()+
  theme(legend.position = "top",   legend.justification='left',
        legend.title = element_blank(),
        plot.title = element_text(color="black", size = 15, face = "bold"),
        axis.title.x = element_text(color="black", size=14, face="bold"),
        axis.text.y = element_text(color= "black", size=11, face="bold"),
        axis.text.x = element_text(color= "black", size=11, face="bold"),
        axis.title.y = element_text(color="black", size=14, face="bold")) +
  expand_limits(y=c(0,14))


plot8 <- plot8 + 
  ggplot2::annotate("text", x = ymd("2020-05-25"), y = 13, fontface=2, size=3, 
                    label = "Killing of George Floyd") +
  ggplot2::annotate("text", x = ymd("2020-03-13"), y = 13, fontface=2, size=3, 
                    label = "Trump declared\n a national emergency") 
plot8 



jpeg("blm.jpeg", units="in", width=10, height=7, res=300)

plot8


dev.off()


######################################################################################################
#####  TABLE E4: Full Seasonality Test
######################################################################################################

## First, re-define the time window. create two different datasets 

## Time Window from 2019 JUNE TO 2020 JUNE #################### 

alldata <- readRDS('ps_alldata.rds')

data_weekly <- alldata
data_weekly$date <- as.Date(data_weekly$date)
data_weekly$date  <- ymd(data_weekly$date)
saturdays <- data_weekly[wday(data_weekly$date) == 7, ] # filter for Saturdays
startDate <- min(saturdays$date) # select first Saturday
data_weekly$week <- floor(as.numeric(difftime(data_weekly$date, startDate, units = "weeks")))

data19 <- data_weekly %>% filter(as.Date(date)>=as.Date('2019-01-01') & as.Date(date)<as.Date('2019-06-01'))
data19$placebocutoff  <-  0
data19$placebocutoff[as.Date(data19$date) >= '2019-03-13'] <- 1 

profw_data19 <- data19 %>% 
  group_by(screen_name,gender,rank,placebocutoff, week) %>% 
  dplyr::summarise(total=n(),family=mean(family,na.rm=T),work=mean(work,na.rm=T))

profw_data19$family <- profw_data19$family*100
profw_data19$work  <- profw_data19$work*100

profw_data19$pre_post <- profw_data19$placebocutoff

assistant19 <- profw_data19 %>% filter(rank==1)
associate19 <- profw_data19 %>% filter(rank==2)
full19 <- profw_data19 %>% filter(rank==3)


data20 <- data_weekly %>% filter(as.Date(date)>=as.Date('2020-01-01') & as.Date(date)<as.Date('2020-06-01'))

profw_data20  <- data20 %>% 
  group_by(screen_name,gender,rank,pre_post, week) %>% 
  dplyr::summarise(total=n(),family=mean(family,na.rm=T),work=mean(work,na.rm=T))

profw_data20$family <- profw_data20$family*100
profw_data20$work <- profw_data20$work*100



assistant20 <- profw_data20 %>% filter(rank==1)
associate20 <- profw_data20 %>% filter(rank==2)
full20 <- profw_data20 %>% filter(rank==3)



######################## Seasonality Test #######################################
##### compare Jan 1 - June 1, 2019 with 3/13, 2019 placebo cutoff with the 2020 data  
##########################################################################################

################### TABLE 2019 with placebo cutoff ###################

m1 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = profw_data19,   exactDOF=TRUE,psdef=FALSE)

m2 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = profw_data19,   exactDOF=TRUE,psdef=FALSE)


m3 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = assistant19,   exactDOF=TRUE,psdef=FALSE)

m4 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data =assistant19,   exactDOF=TRUE,psdef=FALSE)

m5 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = associate19,   exactDOF=TRUE,psdef=FALSE)

m6 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data =associate19,   exactDOF=TRUE,psdef=FALSE)

m7 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = full19,   exactDOF=TRUE,psdef=FALSE)

m8 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data =full19,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(m1,m2,m3,m4,m5,m6, m7, m8,
                                   dep.var.labels=rep(c('Family','Work'),6),
                                   covariate.labels = c('Female*Pandemic'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   column.labels = c('All Faculty','Assistant','Associate', 'Tenured'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   omit = 'week',
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

################### TABLE 2020 with actual cutoff ###################

m1 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = profw_data20,   exactDOF=TRUE,psdef=FALSE)

m2 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = profw_data20,   exactDOF=TRUE,psdef=FALSE)


m3 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = assistant20,   exactDOF=TRUE,psdef=FALSE)

m4 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data =assistant20,   exactDOF=TRUE,psdef=FALSE)

m5 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = associate20,   exactDOF=TRUE,psdef=FALSE)

m6 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data =associate20,   exactDOF=TRUE,psdef=FALSE)

m7 <- felm(family ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data = full20,   exactDOF=TRUE,psdef=FALSE)

m8 <- felm(work ~ gender:pre_post |  factor(screen_name) + factor(week) | 0 | 
             screen_name, data =full20,   exactDOF=TRUE,psdef=FALSE)


table <- capture.output({stargazer(m1,m2,m3,m4,m5,m6, m7, m8,
                                   dep.var.labels=rep(c('Family','Work'),6),
                                   covariate.labels = c('Female*Pandemic'),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   column.labels = c('All Faculty','Assistant','Associate', 'Tenured'), column.separate = c(2,2,2,2),
                                   no.space=TRUE,
                                   omit = 'week',
                                   star.char = c("*", "**", "***"), 
                                   notes = c("* p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=3,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)









